## Read data
Data <- read.table("HMS_Data.raw", header=T)

## Compute observed strata probabilities

# Question 3
obs_Q3 <- ifelse(is.na(Data$CQ3a_Clinic_or_Hospitals)==F & is.na(Data$FQ3a_Clinic)==F, 1, 0)
Q3_total <- subset(Data, obs_Q3 == 1)
Q3_00 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 0 & Q3_total$FQ3a_Clinic ==0, 1, 0)
Q3_01 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 0 & Q3_total$FQ3a_Clinic ==1, 1, 0)
Q3_10 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 1 & Q3_total$FQ3a_Clinic ==0, 1, 0)
Q3_11 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 1 & Q3_total$FQ3a_Clinic ==1, 1, 0)
Pr00_Q3 <- sum(Q3_00)/sum(obs_Q3)
Pr01_Q3 <- sum(Q3_01)/sum(obs_Q3)
Pr10_Q3 <- sum(Q3_10)/sum(obs_Q3)
Pr11_Q3 <- sum(Q3_11)/sum(obs_Q3)
P.3 <- c(Pr00_Q3, Pr01_Q3, Pr10_Q3, Pr11_Q3)

# Question 4a
obs_Q4a <- ifelse(is.na(Data$CQ4a_Ed_Prim_or_Secondary)==F & is.na(Data$FQ4a_Ed_Prim)==F, 1, 0)
Q4a_total <- subset(Data, obs_Q4a == 1)
Q4a_00 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 0 & Q4a_total$FQ4a_Ed_Prim ==0, 1, 0)
Q4a_01 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 0 & Q4a_total$FQ4a_Ed_Prim ==1, 1, 0)
Q4a_10 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 1 & Q4a_total$FQ4a_Ed_Prim ==0, 1, 0)
Q4a_11 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 1 & Q4a_total$FQ4a_Ed_Prim ==1, 1, 0)
Pr00_Q4a <- sum(Q4a_00)/sum(obs_Q4a)
Pr01_Q4a <- sum(Q4a_01)/sum(obs_Q4a)
Pr10_Q4a <- sum(Q4a_10)/sum(obs_Q4a)
Pr11_Q4a <- sum(Q4a_11)/sum(obs_Q4a)
P.4a <- c(Pr00_Q4a, Pr01_Q4a, Pr10_Q4a, Pr11_Q4a)

# Question 7b
obs_Q7b <- ifelse(is.na(Data$CQ7b_Tr_Roads_or_Services)==F & is.na(Data$FQ7b_Tr_Roads)==F, 1, 0)
Q7b_total <- subset(Data, obs_Q7b == 1)
Q7b_00 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 0 & Q7b_total$FQ7b_Tr_Roads ==0, 1, 0)
Q7b_01 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 0 & Q7b_total$FQ7b_Tr_Roads ==1, 1, 0)
Q7b_10 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 1 & Q7b_total$FQ7b_Tr_Roads ==0, 1, 0)
Q7b_11 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 1 & Q7b_total$FQ7b_Tr_Roads ==1, 1, 0)
Pr00_Q7b <- sum(Q7b_00)/sum(obs_Q7b)
Pr01_Q7b <- sum(Q7b_01)/sum(obs_Q7b)
Pr10_Q7b <- sum(Q7b_10)/sum(obs_Q7b)
Pr11_Q7b <- sum(Q7b_11)/sum(obs_Q7b)
P.7b <- c(Pr00_Q7b, Pr01_Q7b, Pr10_Q7b, Pr11_Q7b)

# Question 7c
obs_Q7c <- ifelse(is.na(Data$CQ7c_Tr_Local_or_Links)==F & is.na(Data$FQ7c_Tr_Local)==F, 1, 0)
Q7c_total <- subset(Data, obs_Q7c == 1)
Q7c_00 <- ifelse(Q7c_total$CQ7c_Tr_Local_or_Links == 0 & Q7c_total$FQ7c_Tr_Local ==0, 1, 0)
Q7c_01 <- ifelse(Q7c_total$CQ7c_Tr_Local_or_Links == 0 & Q7c_total$FQ7c_Tr_Local ==1, 1, 0)
Q7c_10 <- ifelse(Q7c_total$CQ7c_Tr_Local_or_Links == 1 & Q7c_total$FQ7c_Tr_Local ==0, 1, 0)
Q7c_11 <- ifelse(Q7c_total$CQ7c_Tr_Local_or_Links == 1 & Q7c_total$FQ7c_Tr_Local ==1, 1, 0)
Pr00_Q7c <- sum(Q7c_00)/sum(obs_Q7c)
Pr01_Q7c <- sum(Q7c_01)/sum(obs_Q7c)
Pr10_Q7c <- sum(Q7c_10)/sum(obs_Q7c)
Pr11_Q7c <- sum(Q7c_11)/sum(obs_Q7c)
P.7c <- c(Pr00_Q7c, Pr01_Q7c, Pr10_Q7c, Pr11_Q7c)

# Question 11a
obs_Q11a <- ifelse(is.na(Data$CQ11a_Discount_Binary)==F & is.na(Data$FQ11a_Discount_Bin)==F, 1, 0)
Q11a_total <- subset(Data, obs_Q11a == 1)
Q11a_00 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 0 & Q11a_total$FQ11a_Discount_Bin ==0, 1, 0)
Q11a_01 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 0 & Q11a_total$FQ11a_Discount_Bin ==1, 1, 0)
Q11a_10 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 1 & Q11a_total$FQ11a_Discount_Bin ==0, 1, 0)
Q11a_11 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 1 & Q11a_total$FQ11a_Discount_Bin ==1, 1, 0)
Pr00_Q11a <- sum(Q11a_00)/sum(obs_Q11a)
Pr01_Q11a <- sum(Q11a_01)/sum(obs_Q11a)
Pr10_Q11a <- sum(Q11a_10)/sum(obs_Q11a)
Pr11_Q11a <- sum(Q11a_11)/sum(obs_Q11a)
P.11a <- c(Pr00_Q11a, Pr01_Q11a, Pr10_Q11a, Pr11_Q11a)

## Obtain upper and lower bounds
## read libraries
library(lpSolve)

## define functions 

# Assumption 1: randomization
# Assumption 3: informative measurement
# Assumption 4: no misclassification for compliant groups
# Assumption 5: no misclassification when leaders' preferences agree with group decisions

# Assumptions 1, 3 and 5
bound.a135 <- function(P,dir,e){
evpoints <- round(1/e) + 1
q <- 0
bound <- rep(NA, evpoints)
pi.cc <- rep(NA, evpoints)
pi.ac <- rep(NA, evpoints)
pi.aa <- rep(NA, evpoints)
pi.nc <- rep(NA, evpoints)
pi.nn <- rep(NA, evpoints)
pi.dc <- rep(NA, evpoints)
pi.da <- rep(NA, evpoints)
pi.dn <- rep(NA, evpoints)
pi.dd <- rep(NA, evpoints)
status <- rep(NA, evpoints)
Q <- rep(NA, evpoints)
Px0 <- P[1] + P[3]
Px1 <- P[2] + P[4]
for(i in 1:evpoints){
    q <- (i-1)*e
    f.obj <- f.obj <- c(1,0,0,0,0,-1,-1,-1,-1)
    f.con <- matrix(c(1-q,0,0,1-q,1,0,0,q,q,
            0,0,0,q,0,q,q,0,0,
            0,1-q,0,0,0,1-q,0,1-q,0,
            q,q,1,0,0,0,1-q,0,1-q,
            0,0,1,0,1,0,1,1,2,
            1,0,0,0,0,0,0,0,0,
            0,1,0,0,0,0,0,0,0,
            0,0,1,0,0,0,0,0,0,
            0,0,0,1,0,0,0,0,0,
            0,0,0,0,1,0,0,0,0,
            0,0,0,0,0,1,0,0,0,
            0,0,0,0,0,0,1,0,0,
            0,0,0,0,0,0,0,1,0,
            0,0,0,0,0,0,0,0,1), nrow=14, byrow=T)
    f.dir <- c("=","=","=","=","<",
        ">=",">=",">=",">=",">=",">=",">=",">=",">=")
    f.rhs <- c(P,1,0,0,0,0,0,0,0,0,0)
    r <- lp(dir, f.obj, f.con, f.dir, f.rhs)
    bound[i] <- r$objval
    pi.cc[i] <- r$solution[1]
    pi.ac[i] <- r$solution[2]
    pi.aa[i] <- r$solution[3]
    pi.nc[i] <- r$solution[4]
    pi.nn[i] <- r$solution[5]
    pi.dc[i] <- r$solution[6]
    pi.da[i] <- r$solution[7]
    pi.dn[i] <- r$solution[8]
    pi.dd[i] <- r$solution[9]
    status[i] <- r$status
    Q[i] <- q
    }
return(as.data.frame(cbind(Q, bound, status, pi.cc, pi.ac, pi.aa,
    pi.nc, pi.nn, pi.dc, pi.da, pi.dn, pi.dd)))
}

# Assumptions 1, 3 and 4
bound.a134 <- function(P,dir,e){
evpoints <- round(1/e) + 1
q <- 0
bound <- rep(NA, evpoints)
pi.cc <- rep(NA, evpoints)
pi.ac <- rep(NA, evpoints)
pi.aa <- rep(NA, evpoints)
pi.an <- rep(NA, evpoints)
pi.ad <- rep(NA, evpoints)
pi.nc <- rep(NA, evpoints)
pi.na <- rep(NA, evpoints)
pi.nn <- rep(NA, evpoints)
pi.nd <- rep(NA, evpoints)
pi.dc <- rep(NA, evpoints)
pi.da <- rep(NA, evpoints)
pi.dn <- rep(NA, evpoints)
pi.dd <- rep(NA, evpoints)
status <- rep(NA, evpoints)
Q <- rep(NA, evpoints)
Px0 <- P[1] + P[3]
Px1 <- P[2] + P[4]
for(i in 1:evpoints){
    q <- (i-1)*e
    f.obj <- c(1,0,0,0,0,0,0,0,0,-1,-1,-1,-1)
    f.con <- matrix(c(1-q,0,0,0,0,1-q,0,1,q,0,0,q,q,
            0,0,0,0,0,q,1,0,1-q,q,q,0,0,
            0,1-q,0,1,q,0,0,0,0,1-q,0,1-q,0,
            q,q,1,0,1-q,0,0,0,0,0,1-q,0,1-q,
            0,0,1,1,2,0,1,1,2,0,1,1,2,
            1,0,0,0,0,0,0,0,0,0,0,0,0,
            0,1,0,0,0,0,0,0,0,0,0,0,0,
            0,0,1,0,0,0,0,0,0,0,0,0,0,
            0,0,0,1,0,0,0,0,0,0,0,0,0,
            0,0,0,0,1,0,0,0,0,0,0,0,0,
            0,0,0,0,0,1,0,0,0,0,0,0,0,
            0,0,0,0,0,0,1,0,0,0,0,0,0,
            0,0,0,0,0,0,0,1,0,0,0,0,0,
            0,0,0,0,0,0,0,0,1,0,0,0,0,
            0,0,0,0,0,0,0,0,0,1,0,0,0,
            0,0,0,0,0,0,0,0,0,0,1,0,0,
            0,0,0,0,0,0,0,0,0,0,0,1,0,
            0,0,0,0,0,0,0,0,0,0,0,0,1), nrow=18, byrow=T)
    f.dir <- c("=","=","=","=","<",
        ">=",">=",">=",">=",">=",
        ">=",">=",">=",">=",">=",">=",">=",">=")
    f.rhs <- c(P,1,0,0,0,0,0,0,0,0,0,0,0,0,0)
    r <- lp(dir, f.obj, f.con, f.dir, f.rhs)
    bound[i] <- r$objval
    pi.cc[i] <- r$solution[1]
    pi.ac[i] <- r$solution[2]
    pi.aa[i] <- r$solution[3]
    pi.an[i] <- r$solution[4]
    pi.ad[i] <- r$solution[5]
    pi.nc[i] <- r$solution[6]
    pi.na[i] <- r$solution[7]
    pi.nn[i] <- r$solution[8]
    pi.nd[i] <- r$solution[9]
    pi.dc[i] <- r$solution[10]
    pi.da[i] <- r$solution[11]
    pi.dn[i] <- r$solution[12]
    pi.dd[i] <- r$solution[13]
    status[i] <- r$status
    Q[i] <- q
    }
return(as.data.frame(cbind(Q, bound, status, pi.cc,
    pi.ac, pi.aa, pi.an, pi.ad, pi.nc, pi.na, pi.nn, pi.nd,
    pi.dc, pi.da, pi.dn, pi.dd)))
}


# Assumptions 1 and 3
bound.a13 <- function(P,dir,e){
evpoints <- round(1/e) + 1
q <- 0
bound <- rep(NA, evpoints)
pi.cc <- rep(NA, evpoints)
pi.ca <- rep(NA, evpoints)
pi.cn <- rep(NA, evpoints)
pi.cd <- rep(NA, evpoints)
pi.ac <- rep(NA, evpoints)
pi.aa <- rep(NA, evpoints)
pi.an <- rep(NA, evpoints)
pi.ad <- rep(NA, evpoints)
pi.nc <- rep(NA, evpoints)
pi.na <- rep(NA, evpoints)
pi.nn <- rep(NA, evpoints)
pi.nd <- rep(NA, evpoints)
pi.dc <- rep(NA, evpoints)
pi.da <- rep(NA, evpoints)
pi.dn <- rep(NA, evpoints)
pi.dd <- rep(NA, evpoints)
status <- rep(NA, evpoints)
Q <- rep(NA, evpoints)
Px0 <- P[1] + P[3]
Px1 <- P[2] + P[4]
for(i in 1:evpoints){
    q <- (i-1)*e
    f.obj <- c(1,1,1,1,0,0,0,0,0,0,0,0,-1,-1,-1,-1)
    f.con <- matrix(c(1-q,0,1-q,0,0,0,0,0,1-q,0,1,q,0,0,q,q,
            0,1-q,0,1-q,0,0,0,0,q,1,0,1-q,q,q,0,0,
            0,0,q,q,1-q,0,1,q,0,0,0,0,1-q,0,1-q,0,
            q,q,0,0,q,1,0,1-q,0,0,0,0,0,1-q,0,1-q,
            0,1,1,2,0,1,1,2,0,1,1,2,0,1,1,2,
            1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
            0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
            0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
            0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
            0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
            0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,
            0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
            0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
            0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
            0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
            0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
            0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
            0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
            0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
            0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1), nrow=21, byrow=T)
    f.dir <- c("=","=","=","=","<",
        ">=",">=",">=",">=",">=",">=",">=",">=",
        ">=",">=",">=",">=",">=",">=",">=",">=")
    f.rhs <- c(P,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
    r <- lp(dir, f.obj, f.con, f.dir, f.rhs)
    bound[i] <- r$objval
    pi.cc[i] <- r$solution[1]
    pi.ca[i] <- r$solution[2]
    pi.cn[i] <- r$solution[3]
    pi.cd[i] <- r$solution[4]
    pi.ac[i] <- r$solution[5]
    pi.aa[i] <- r$solution[6]
    pi.an[i] <- r$solution[7]
    pi.ad[i] <- r$solution[8]
    pi.nc[i] <- r$solution[9]
    pi.na[i] <- r$solution[10]
    pi.nn[i] <- r$solution[11]
    pi.nd[i] <- r$solution[12]
    pi.dc[i] <- r$solution[13]
    pi.da[i] <- r$solution[14]
    pi.dn[i] <- r$solution[15]
    pi.dd[i] <- r$solution[16]
    status[i] <- r$status
    Q[i] <- q
    }
return(as.data.frame(cbind(Q, bound, status, pi.cc, pi.ca, pi.cn, pi.cd,
    pi.ac, pi.aa, pi.an, pi.ad, pi.nc, pi.na, pi.nn, pi.nd,
    pi.dc, pi.da, pi.dn, pi.dd)))
}

by <- 0.0005

# Question 3
resQ3.135.u <- bound.a135(P.3, "max", by)
resQ3.135.l <- bound.a135(P.3, "min", by)
resQ3.134.u <- bound.a134(P.3, "max", by)
resQ3.134.l <- bound.a134(P.3, "min", by)
resQ3.13.u <- bound.a13(P.3, "max", by)
resQ3.13.l <- bound.a13(P.3, "min", by)

# Question 4a
resQ4a.135.u <- bound.a135(P.4a, "max", by)
resQ4a.135.l <- bound.a135(P.4a, "min", by)
resQ4a.134.u <- bound.a134(P.4a, "max", by)
resQ4a.134.l <- bound.a134(P.4a, "min", by)
resQ4a.13.u <- bound.a13(P.4a, "max", by)
resQ4a.13.l <- bound.a13(P.4a, "min", by)

# Question 7b
resQ7b.135.u <- bound.a135(P.7b, "max", by)
resQ7b.135.l <- bound.a135(P.7b, "min", by)
resQ7b.134.u <- bound.a134(P.7b, "max", by)
resQ7b.134.l <- bound.a134(P.7b, "min", by)
resQ7b.13.u <- bound.a13(P.7b, "max", by)
resQ7b.13.l <- bound.a13(P.7b, "min", by)

# Question 7c
resQ7c.135.u <- bound.a135(P.7c, "max", by)
resQ7c.135.l <- bound.a135(P.7c, "min", by)
resQ7c.134.u <- bound.a134(P.7c, "max", by)
resQ7c.134.l <- bound.a134(P.7c, "min", by)
resQ7c.13.u <- bound.a13(P.7c, "max", by)
resQ7c.13.l <- bound.a13(P.7c, "min", by)

# Question 11a
resQ11a.135.u <- bound.a135(P.11a, "max", by)
resQ11a.135.l <- bound.a135(P.11a, "min", by)
resQ11a.134.u <- bound.a134(P.11a, "max", by)
resQ11a.134.l <- bound.a134(P.11a, "min", by)
resQ11a.13.u <- bound.a13(P.11a, "max", by)
resQ11a.13.l <- bound.a13(P.11a, "min", by)

## Exclude lines where bounds are not identified
Q <- seq(0, 1, by=by)

Q3.135.ubound <- ifelse(resQ3.135.u$bound > resQ3.135.l$bound, resQ3.135.u$bound, NA)
Q3.135.lbound <- ifelse(resQ3.135.u$bound > resQ3.135.l$bound, resQ3.135.l$bound, NA)
Q4a.135.ubound <- ifelse(resQ4a.135.u$bound > resQ4a.135.l$bound, resQ4a.135.u$bound, NA)
Q4a.135.lbound <- ifelse(resQ4a.135.u$bound > resQ4a.135.l$bound, resQ4a.135.l$bound, NA)
Q7b.135.ubound <- ifelse(resQ7b.135.u$bound > resQ7b.135.l$bound, resQ7b.135.u$bound, NA)
Q7b.135.lbound <- ifelse(resQ7b.135.u$bound > resQ7b.135.l$bound, resQ7b.135.l$bound, NA)
Q7c.135.ubound <- ifelse(resQ7c.135.u$bound > resQ7c.135.l$bound, resQ7c.135.u$bound, NA)
Q7c.135.lbound <- ifelse(resQ7c.135.u$bound > resQ7c.135.l$bound, resQ7c.135.l$bound, NA)
Q11a.135.ubound <- ifelse(resQ11a.135.u$bound > resQ11a.135.l$bound, resQ11a.135.u$bound, NA)
Q11a.135.lbound <- ifelse(resQ11a.135.u$bound > resQ11a.135.l$bound, resQ11a.135.l$bound, NA)

Q3.134.ubound <- ifelse(resQ3.134.u$bound > resQ3.134.l$bound, resQ3.134.u$bound, NA)
Q3.134.lbound <- ifelse(resQ3.134.u$bound > resQ3.134.l$bound, resQ3.134.l$bound, NA)
Q4a.134.ubound <- ifelse(resQ4a.134.u$bound > resQ4a.134.l$bound, resQ4a.134.u$bound, NA)
Q4a.134.lbound <- ifelse(resQ4a.134.u$bound > resQ4a.134.l$bound, resQ4a.134.l$bound, NA)
Q7b.134.ubound <- ifelse(resQ7b.134.u$bound > resQ7b.134.l$bound, resQ7b.134.u$bound, NA)
Q7b.134.lbound <- ifelse(resQ7b.134.u$bound > resQ7b.134.l$bound, resQ7b.134.l$bound, NA)
Q7c.134.ubound <- ifelse(resQ7c.134.u$bound > resQ7c.134.l$bound, resQ7c.134.u$bound, NA)
Q7c.134.lbound <- ifelse(resQ7c.134.u$bound > resQ7c.134.l$bound, resQ7c.134.l$bound, NA)
Q11a.134.ubound <- ifelse(resQ11a.134.u$bound > resQ11a.134.l$bound, resQ11a.134.u$bound, NA)
Q11a.134.lbound <- ifelse(resQ11a.134.u$bound > resQ11a.134.l$bound, resQ11a.134.l$bound, NA)

Q3.13.ubound <- ifelse(resQ3.13.u$bound > resQ3.13.l$bound, resQ3.13.u$bound, NA)
Q3.13.lbound <- ifelse(resQ3.13.u$bound > resQ3.13.l$bound, resQ3.13.l$bound, NA)
Q4a.13.ubound <- ifelse(resQ4a.13.u$bound > resQ4a.13.l$bound, resQ4a.13.u$bound, NA)
Q4a.13.lbound <- ifelse(resQ4a.13.u$bound > resQ4a.13.l$bound, resQ4a.13.l$bound, NA)
Q7b.13.ubound <- ifelse(resQ7b.13.u$bound > resQ7b.13.l$bound, resQ7b.13.u$bound, NA)
Q7b.13.lbound <- ifelse(resQ7b.13.u$bound > resQ7b.13.l$bound, resQ7b.13.l$bound, NA)
Q7c.13.ubound <- ifelse(resQ7c.13.u$bound > resQ7c.13.l$bound, resQ7c.13.u$bound, NA)
Q7c.13.lbound <- ifelse(resQ7c.13.u$bound > resQ7c.13.l$bound, resQ7c.13.l$bound, NA)
Q11a.13.ubound <- ifelse(resQ11a.13.u$bound > resQ11a.13.l$bound, resQ11a.13.u$bound, NA)
Q11a.13.lbound <- ifelse(resQ11a.13.u$bound > resQ11a.13.l$bound, resQ11a.13.l$bound, NA)

## Create plots

pdf("fig1.pdf", width = 11, height = 8.5)

par(cex=1.45)
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q7c: Villages (0) or Major Centers (1)?",
    xlab = expression("Treatment Assignment Probability"~ ~Q == Pr(Z[i]^"*" == 1)),
    ylab = "Average Treatment Effect")
lines(Q, Q7c.13.ubound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q7c.13.lbound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q7c.134.ubound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q7c.134.lbound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q7c.135.ubound, col="darkgreen", lwd=2)
lines(Q, Q7c.135.lbound, col="darkgreen", lwd=2)
abline(0,0)
points(P.7c[2]+P.7c[4], 
    P.7c[4]/(P.7c[2]+P.7c[4]) - P.7c[3]/(P.7c[1]+P.7c[3]),
    pch = 19)
legend(0, 1.15, c("Persuasion Scenario", "Incentive Scenario",
    "No Substantive Assumption"), box.lty = 0, lwd = c(2,3,4),
    lty = c(1,2,3), col = c("darkgreen", "darkred", "darkblue"))
text(0.30, 0.15, expression(frac(P[0*0],1-Q) - frac(P[0*1],Q)), cex=0.9)
text(0.65, 0.15, expression(frac(P[1*1],Q) - frac(P[1*0],1-Q)), cex=0.9)
text(0.1, -0.95, expression(- frac(P[1*0]+P[1*1],1-Q)), cex=0.9)
text(0.4, -0.5, expression(- frac(P[0*1],Q) - frac(P[1*0],1-Q)), cex=0.9)
text(0.8, -0.85, expression(- frac(P[0*0]+P[0*1],Q)), cex=0.9)
arrows(0.375, 0.15, 0.425, 0.0875, length = 0.1)
arrows(0.575, 0.15, 0.535, 0.089, length = 0.1)
arrows(0.175, -0.95, 0.205, -0.775, length = 0.1)
arrows(0.4, -0.6, 0.39, -0.78, length = 0.1)
arrows(0.735, -0.85, 0.67, -0.6, length = 0.1)

dev.off()


pdf("fig2.pdf", width = 11, height = 8.5)
par(mfrow=c(2,2))

# Question 3
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q3: Clinics (0) or Hospitals (1)?",
    xlab = expression("Treatment Assignment Probability"~ ~Q == Pr(Z[i]^"*" == 1)),
    ylab = "Average Treatment Effect")
lines(Q, Q3.13.ubound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q3.13.lbound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q3.134.ubound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q3.134.lbound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q3.135.ubound, col="darkgreen", lwd=2)
lines(Q, Q3.135.lbound, col="darkgreen", lwd=2)
abline(0,0)
points(P.3[2]+P.3[4], 
    P.3[4]/(P.3[2]+P.3[4]) - P.3[3]/(P.3[1]+P.3[3]),
    pch = 19)
    
# Question 4a
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q4a: Primary (0) or Secondary (1) Education?",
    xlab = expression("Treatment Assignment Probability"~ ~Q == Pr(Z[i]^"*" == 1)),
    ylab = "Average Treatment Effect")
lines(Q, Q4a.13.ubound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q4a.13.lbound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q4a.134.ubound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q4a.134.lbound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q4a.135.ubound, col="darkgreen", lwd=2)
lines(Q, Q4a.135.lbound, col="darkgreen", lwd=2)
abline(0,0)
points(P.4a[2]+P.4a[4], 
    P.4a[4]/(P.4a[2]+P.4a[4]) - P.4a[3]/(P.4a[1]+P.4a[3]),
    pch = 19)

# Question 7b
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q7b: Roads (0) or Public Transportation (1)?",
    xlab = expression("Treatment Assignment Probability"~ ~Q == Pr(Z[i]^"*" == 1)),
    ylab = "Average Treatment Effect")
lines(Q, Q7b.13.ubound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q7b.13.lbound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q7b.134.ubound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q7b.134.lbound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q7b.135.ubound, col="darkgreen", lwd=2)
lines(Q, Q7b.135.lbound, col="darkgreen", lwd=2)
abline(0,0)
points(P.7b[2]+P.7b[4], 
    P.7b[4]/(P.7b[2]+P.7b[4]) - P.7b[3]/(P.7b[1]+P.7b[3]),
    pch = 19)

# Question 11a
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q11a: Consume (0) or Invest (1) Windfalls?",
    xlab = expression("Treatment Assignment Probability"~ ~Q == Pr(Z[i]^"*" == 1)),
    ylab = "Average Treatment Effect")
lines(Q, Q11a.13.ubound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q11a.13.lbound, col="darkblue", lty="dotted", lwd=4)
lines(Q, Q11a.134.ubound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q11a.134.lbound, col="darkred", lty="dashed", lwd=3)
lines(Q, Q11a.135.ubound, col="darkgreen", lwd=2)
lines(Q, Q11a.135.lbound, col="darkgreen", lwd=2)
abline(0,0)
points(P.11a[2]+P.11a[4], 
    P.11a[4]/(P.11a[2]+P.11a[4]) - P.11a[3]/(P.11a[1]+P.11a[3]),
    pch = 19)

dev.off()

